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£Sj . Abstract 

The intrinsically nonlinear nature of quantum field theory provides a fundamental 
■ complication for lattice calculations, when the physical implications of the subtleties 

of Fourier theory are taken into account. Even though the fundamental fields are con- 
| strained to the first Brillouin zone, Fourier theory tells us that the high-momentum 

components of products of these fields "bleed into" neighbouring Brillouin zones, where 
they "alias" (or "masquerade") as low-momentum contributions, violating the conser- 
vation of energy and momentum, and fundamentally distorting calculations. In this 
paper I offer a general strategy for eliminating the artefacts of aliasing in practical 
calculations. 
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1. Introduction and motivation 



There is a strange relationship between the ways that engineers and theoretical physicists 
employ Fourier theory. Engineers generally use it for physical applications that are, from the 
theoretical physicist's point of view, almost trivally simple in structure: physical equations 
that are often linear, or, at worst, nonlinear in ways that are simple to understand, and 
even simpler to describe mathematically. Physicists, in contrast, analyse almost intractably 
complicated mathematical descriptions of physical reality, for which Fourier techniques are 
but one of the fundamental tools in what can be a vast mathematical toolbox of almost 
unbelievable sophistication and abstraction. 

It is ironic, then, that the average engineer often gets a more thorough grounding in the 
fundamentals of Fourier theory than the average theoretical physicist. I don't know why 
this is so; it seems to be something of a historical accident. That the same mathematical 
formalism has been developed, independently, in both fields is illustrated most starkly by 
the fact that the fundamental theorems and constructs of Fourier theory are named after 
different people, depending on which faculty department is teaching it! 

It is perhaps a truism that engineers can spend much more time analysing every sub- 
tlety of the physical theory they are using to describe their creations, simply because such 
descriptions are so fundamentally simple and straightforward (from the physicist's point 
of view, anyway). Nevertheless, it seems to me that there are some lessons that engineers 
have long learnt, from their extensive application of Fourier theory to real- world applications 
(when a simple mathematical oversight can mean the difference between a device working as 



1 



designed or failing dismally), that have not been sufficiently hammered home to theoretical 
physicists — if they have, indeed, been explicitly recognised as problems at all. 

1 believe that the most insidious of these is rooted in the simple process of forming 
products of fields in lattice field theory. It is a fundamental theorem of Fourier theory 
that forming the product of two fields in position space leads to a "convolution" of their 
momentum-space representations. This process axiomatically leads to the Fourier trans- 
form of the product "bleeding out" of what the physicist calls the "first Brillouin zone" 
(the engineer calls it "above the Nyquist frequency") into the surrounding "zones", which 
necessarily means that it is automatically "shifted down" in momentum. The engineer calls 
this phenomenon "aliasing" : high-momentum components "masquerade" as low-momentum 
components, and generally completely destroy the fidelity of the low-frequency signal in the 
process. The physicist generally describes the result as being analogous to an "Umklapp 
process" in a crystal, because this is the physical example (with real lattices, no less) for 
which this phenomenon is most familiar. 

No matter what it is called, this phenomenon can be insidiously devastating for any 
lattice calculation aiming to obtain as accurate a result as possible for a given amount of 
computing power. The fact that this phenomenon violates the conservation of energy and 
momentum should be sufficient to ring alarm bells for any theorist (there is no physically real 
"crystal" to supply the Umklapp momentum in quantum field theory — it just comes from 
nowhere!). That fields end up being in the wrong place in momentum space is of concern 
not just to the theorist, but also the pragmatist who simply wants to extract physically 
meaningful numbers from a lattice calculation. 

A simple example will suffice to demonstrate the general havoc wreaked by aliasing. 
Imagine that we wish to compute some sort of expectation value, that depends on a product 
of some number of fundamental fields, together with a number of operators acting on the 
fields. Computing the expectation value corresponds to evaluating the Fourier transform 
of this result at zero momentum. However, if aliasing is not prevented, we find that there 
are two contributions to the result. The first is the true zero- momentum result due to the 
latticised approximation of the physical system (with whatever unavoidable approximations 
and inaccuracies that that entails), which is what we are trying to extract. The second is 
due to the interaction of the components of each field (and the operators in question) at high 
momentum values, which has been aliased down to zero momentum due to the fundamental 
properties of Fourier theory. The result, of course, comes to us as a single number, with 
these two contributions inextricably intertwined. 

It may come as a shock, to some, that this second contribution is present at all, and is to a 
greater or lesser extent confounding the entire lattice calculation by its presence. Indeed, the 
components that make up the second contribution rightly belong at high momentum values— 
outside the first Brillouin zone completely. What are they doing back at zero momentum? 

Such is the fun and games that one encounters if one entrusts Fourier theory with our 
calculations, without first "reading the fine print" of that mathematical formalism. In the 
rest of this paper, I will try to provide a "user's guide" for Fourier theory, in which this "fine 
print" is enlarged, and made intelligible. For many readers, much of what I will review here 
is like "kindergarten physics" — not worried about since undergraduate days; but when it is 
all put together, the importance of its ramifications for lattice calculations may be surprising. 
Other readers may find useful the collecting together of the fundamental theorems of Fourier 
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theory, and its connection with physical reality, which is often developed somewhat flippantly 
and "on the run" in elementary physics applications. Finally, 1 will offer a general strategy for 
removing the aliasing problem at its root, that may be of use in practical lattice calculations. 



2. A review of the fundamentals of Fourier theory 



In this section 1 provide a brief overview of the fundamental mathematical elements of 
Fourier theory. 1 concentrate on one-dimensional lattices; the extension to an arbitrary num- 
ber of dimensions is elementary, and does not introduce any new constructs. The physical 
application and interpretation of this purely mathematical machinery will be discussed in 
the next section. 

We start with a set of iV lattice sites. We can label these sites with an integral index n 
that falls within a range of iV successive integers. For example, we might have n running 
from to N— 1; but any other such range is equally possible. For definiteness, for the rest 
of this section we will assume that n runs from 1 to N. 

We can define a function, say, f n , that is defined on each of the N lattice sites. The N 
values f n are not restricted in any way at all. They may follow some mathematical rule— 
and in practical applications usually will — but in principle they may be simply N numbers 
created at random. 

We now define the Fourier transform, F b = J-'ifn}, as follows: 



The values F b represent lattice sites in "Fourier space". The index b, like n, can fall within 
any range of N successive integers, which need not match the range of n. For definiteness, 
however, for the rest of this section we will assume that b also runs from 1 to N. 

The fundamental theorem of Fourier theory is that the function f n can be retrieved from 
the function F b by means of the inverse Fourier transform, f n = 



which differs from (1) only in that n and b have been interchanged, / and F have been 
interchanged, and —% has been replaced by +i in the complex exponential. It is instructive 
to review the proof of the fundamental theorem, which is straightforward. We do this by 
evaluating the right-hand side of (2). Substituting F b from (1), we have 




2-Kibn/N r 
Jn- 





b- 



(2) 




2wibn'/N j ; 




N 



1 



(3) 
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where in the last line we have simply interchanged the order of performing the sums, which 
is permissible because the sums are manifestly finite. We now consider the last sum in (3), 



N 

2irib(n-n')/N 



e 



E 

6=1 

There are two possible cases. If n = n', the exponent will vanish for any value of b, and hence 
the complex exponential is unity for all values of b, and so we are simply summing up N 
copies of unity, which is just N. If n ^ n', on the other hand, it is straightforward to verify 
that the sum vanishes by symmetry. In the case that \n — n'\ is coprime to N, then we are 
simply adding together all the N-th roots of unity, which are equally spaced around the unit 
circle in the Argand plane, and hence sum (vectorially) to zero. In the alternative case that 
\n — n'\ is not coprime to N, then if g — gcd(|n — n'\, N), we are simply adding together all 
the (N/g)-th roots of unity, a total of g times. Thus 

N 

Y,e 2mb{n - n ' )/N =N5 nn ,, (4) 

6=1 

where is the Kronecker delta symbol. Substituting (4) into (3) yields 

1 N 

J-j: E fn'NSnn' = f n , 
n'=l 

hence establishing the fundamental theorem. The theorem equally provides that the F\, may 
be uniquely recovered from the /„. Indeed, there is no fundamental mathematical distinction 
between the "original" space and "Fourier" space; they simply provide two distinct but 
equivalent representations of iV complex quantities. 

The expressions for the Fourier transformation (1) and its inverse (2) show why the 
indices n and b may be shifted by any arbitrary integral multiple of iV without changing the 
mathematical results. Namely, if we replace n by n + cN, where c is an integer, then the 
complex exponential factor in (2) or (1) becomes 

e ±2nib(n+cN)/N _ e ±2nibn / N e ±2wibc _ e ±2wibn/N 

because e ±2mbc = 1 if b and c are integers, as they are here. Likewise, if we replace b by 
b + cN, then the complex exponential factor in (2) or (1) becomes 

e ±2ni(b+cN)n/N _ e ±2iribn/N e ±2nicn _ e ±2nibn/N 

Let us now consider what happens if we form the product of two functions /„ and g n , 

(fg)n = fn9n- 

If we take the Fourier transform of (fg) n , we have, from (1), 

1 N 

HUgU = -j= E e~ 2mbn/N f n g n . 
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If we now use (2) to write f n and g n in terms of their Fourier transforms Fb and Gb, this 
becomes 

-i N -i N I N 

1 -2iribn/N 1 +2mb'n/N rp 1 +2nib"n/N pi 

N N N 

E E E e 2m{b ' +6 "- 6)n/7V F 5 ,GV 



J v n=l fe'=l b"=l 

1 AT AT AT 

= ^E^E ^" E e W- b W", (5) 

JV b'=l b"=l n=l 

where again we can interchange the order of the sums because they are all manifestly finite. 
We can now use the result (4) for the final sum in (5), to obtain 

^ n n I N 

TT^h E F b' E Gb»N5v>b-v = -7= E F b'G b ~b>- 

iv b'=l b"=l ViV 6 / =1 

This last bilinear operation in F b and Gb is of such fundamental importance that it is given a 
special name, the convolution of F and G, and is denoted by a special symbol, the asterisk: 

1 N 

(F*G) b = ^=J2 F b'Gb-b'. (6) 
VN b , =1 

We have thus shown that the Fourier transform of a product of functions is the convolution 
of their Fourier transforms: 

HifgU = (F*G) b . (7) 

It is straightforward to verify that the same theorem holds true for the inverse Fourier 
transform of a product of functions in Fourier space: 

J 7 ~ 1 {(FG)b} = (f*g) n , (8) 

where convolution in the "original" space takes the exact same form as (6): 

1 N 



(f*g)n = -7^ E fn'9n-n>- 
V N n , =1 



(9) 



Now, the reader with a careful eye may have noticed that there is a slight problem with 
the definitions (6) and (9) of the convolution operation. We have specified, for definiteness, 
that the indices n and b will each take on values from 1 to N, for the purposes of this section. 
The indices in (6) and (9), however, take on values that clearly lie outside this range. For 
example, consider the computation of the 6 = 2 value in the convolution (6): 

1 N 

v N h >=i 

= --j^^FiGi + F 2 G + F 3 G_i + . . . + F N _iG- N+3 + F N G^ N+2 Y 



The index on G should run from 1 to N, like that on F; but here it is running from —N+ 2 
to 1. So what happened? 

The answer is that we made a rather uncritical use of the result (4) in deriving the 
convolution results (7) and (8). The correct generalisation of (4) is actually the following: 



where c is an integer. In the case of the identity (4), we had c = n — n', so that the condition 
that selects the nonzero term is n — n = (mod N), namely n = n' (mod N). But since 
both n and n' took on values between 1 and N, the only way that this congruence could 
hold was for the case of n — n'; hence the result (4) was, in fact, correct. For the final sum 
in (5), on the other hand, the correct application of (4') yields 



In other words, we should not have selected out the term for which b" = b — b' (which may, 
in general, fall outside the range of values taken on by indices b), but rather we should 
have looked for that (permissible) value of b" that was congruent to b — b' modulo N. For 
some values of b and b', the permissible value of b" will simply be b — b'; for others, it will 
be N + b — b'. If we had chosen a different range of values for b, then other multiples of N 
would have to be added in order to obtain values of b" within its allowed range. 

Catering for these various cases in the mathematical definition of the convolution is a 
severe inconvenience, and spoils the elegance of its expression. The usual alternative is to 
deem that all indices n and b are, in fact, permitted to take on any integral values, but 
that it is implicitly assumed that any index arithmetic is ultimately carried out modulo N, 
and then shifted into the chosen range of values for indices of that type. This effectively 
stipulates that the lattices in n-space and 6-space are to be understood to be formed into 
rings of N sites each. 

Whilst this prescription provides a mathematical solution, its physical ramifications are 
both subtle and potentially devastating. Indeed, it is this "trick" that is fundamentally 
responsible for largely hiding the aliasing problem that has affected the fidelity of lattice 
calculations in the past. 



The previous section has reviewed the fundamental elements of Fourier theory in purely 
mathematical terms, without reference to any physical application or interpretation of the 
results. We shall now review how correspondences between these mathematical constructs 
and physical reality are generally constructed. 

The lattice of sites in the "original" n-space is usually taken to correspond to equispaced 
positions in position space; namely, we deem that the position x n corresponding to the lattice 
site n is just 




N if c = (mod iV), 
ifc^O (mod iV), 



(4') 




N if b" = b — b' (mod N), 
if b"^ b -b' (mod N). 



3. Employing lattices for physical applications 



(10) 
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where a is the "lattice spacing", namely, the real (dimensionful) distance between adjacent 
lattice sites. (Recall that we are restricting our attention here to one-dimensional lattices; 
the generalisation to an arbitrary number of dimensions is elementary.) The lattice has a 
total (real) length of 

L = Na. (11) 

The relation (2) now tells us that any arbitrary function f n defined on the n sites x n can 
be written as a suitable linear superposition of iV complex exponential "waves". Based on 
our experience with quantum mechanics, we interpret each such wave as being an "eigenstate 
of momentum". Now, inverting (10) for n, we have 

2-nibn/N _ 2iribx n /Na 



For this to be equivalent to the oscillating part of our usual definition of a plane wave 
eigenstate of momentum p, namely, e tpx (where I shall always use units in which h — 1), we 
therefore require 

» - Wa b - (12 » 

In other words, the momentum value p is simply proportional to the index b. Now, since 
the "length" of the lattice in 6-space is also N sites, the relationship (12) tells us that the 
"length" of the lattice in momentum space is just 

A = -. (13) 
a 

(The somewhat annoying factors of 2ti in (12) and (13), that float around all Fourier ex- 
pressions involving p, fundamentally arise because we have chosen our units such that h — 1 
rather than h — 1. Engineers essentially make the alternative choice, so that their "momen- 
tum space" variable is not complicated in this way. However, this throws an extra factor 
of 27r into the Fourier transform of the derivative operator, and so is inconvenient when we 
want to deal extensively with the formulation of differential equations, as we do in quantum 
field theory.) 

So far, our connection of n-space to position space, and 6-space to momentum space, has 
retained the fundamental mathematical symmetry between the two spaces that was evident 
in Sec. 2. However, once we start thinking about these spaces with a physical interpretation in 
mind, this mathematical symmetry begins to dissolve, because position space and momentum 
space do not enter into the laws of physics in identical ways — far from it, in fact. 

Consider the situation in position space. Our fundamental laws of physics are transla- 
tionally invariant, which means that no particular position x is fundamentally different than 
any other particular position x'. This translational invariance is not obeyed by the position- 
space lattice, as we have defined it so far: it has a finite length L, and hence has two "ends" 
which have only one neighbouring lattice site each. In this context, we willingly embrace the 
"ring" interpretation of the index n described in the previous section, so that position space 
is effectively wrapped into a ring of N lattice sites. This ensures that there are no "ends", 
albeit at the expense of imposing a fundamental periodicity on position space — namely, the 
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variable x can be defined to run from x = — oo to x — +00, but every function of x is forced 
to repeat identically with a periodicity of L. Our gut physical instinct tells us that, if we can 
make L sufficiently large for the physical system we are wishing to model, the "finite volume" 
artefacts of this periodicity (i.e., the interactions between one "copy" of the physical system 
in position space and the adjacent "copies" of itself) can be made sufficiently small that they 
do not detract substantially from the fidelity of the mathematical model we are employing. 

The situation in momentum space is completely different. Here we again have a lattice 
of finite "length" A. However, we do not have any fundamental property of "momentum 
translation invariance" that we wish to maintain. (The relativity of motion of Galilean or 
Lorentzian kinematics is somewhat more subtle than a simple translation in momentum 
space, and we can generally perform a transformation to the "centre of momentum" frame 
of reference, so that such invariance is of no practical concern anyway.) Rather, we con- 
sider the momentum "length" A to provide an upper limit to the momentum states that 
can contribute to our calculations. Indeed, since we have (in any sensible calculation) made 
the abovementioned transformation to a suitable "centre of momentum" frame of reference, 
we make use of our freedom to choose the range of b (and hence p) to ensure that momen- 
tum space is "centred" on p = 0; in other words, we choose to have p falling in the range 
—A < p < +A. In the context of quantum field theory, this "cutoff" in momentum space 
is actually quite convenient, regulating the formalism automatically by removing all of the 
higher-momentum states from our calculations. 

Again, our very language here demonstrates the fundamentally different way that we 
view momentum space compared to position space. We accept implicitly the discussion 
of a "momentum state", because it is something that we are used to considering, even in 
undergraduate physics. Mathematically, however, a "momentum state" is actually quite a 
singular object. It corresponds to a function that has but a single component in momen- 
tum space. For example, if we are considering the momentum state corresponding to the 
momentum-space index 6 = 3, then the function in momentum space takes the form 

F b = A5 b3 , (14) 

where A is some arbitrary amplitude. The reason that we accept such a function, without 
need for explicit comment, is that we have been long accustomed to the fact that any 
translationally-invariant wave equation will have eigenstates that are of precisely this form. 
Indeed, an engineer would equally accept a plane wave function without argument — perhaps, 
at worst, insisting on a real wave, such as 

F b = ^{6 b3 + 5 b ,_ 3 }. (15) 

In contrast, we wouldn't generally talk about a "position state"; such a phrase would need to 
be explained in order to be intelligible. However, if we were to write down the position-space 
analogue of (14), namely, 

F n = A5 n3 , 

then we would immediately recognise what is being considered: a delta-function in position 
space. Why the different description? Simply because such a function is not an eigenstate 
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of any naturally-occurring equation of physics that we are used to considering. Rather, it is 
a fairly advanced mathematical construct that is used in more abstract analyses of physical 
systems. The engineer likewise accepts it as such; but in this case there is no objection to its 
form — no insistence that it be "symmetrised" into a form similar to (15). Why not? Simply 
because it is in position space — "real" space — that the engineer insists on quantities being 
real (and indeed the physicist concurs, with the refinement that only physically observable 
quantities need actually be real in position space). Momentum space, in contrast, is an 
"artificial" or "mathematical" space, with which we learn to make an intellectual connection, 
but which is not the "real" world of our physical senses. Indeed, for a function to be real in 
position space, its Fourier transform needs to be Hermitian (not real) in momentum space, 
namely, 

F. h = F b *, (16) 

where the index —b is understood to be taken modulo N and shifted into the correct range for 
momentum indices (although if we have centred the momentum index range on p = 0, then 
this caveat is not required). If we were to furthermore insist that the function in momentum 
space be real, then we would be restricting the function in position space to be an even 
function of x, which is not appropriate or applicable in general. Thus, the physicist and the 
engineer both accept, without question, that functions in momentum space will, in general, 
be complex-valued, even when they are describing real functions in position space. 

Now, let us leave the divergence between our physical interpretation of position space and 
momentum space to one side, for the moment. Instead, let us focus on a more fundamental 
question: how do we connect the continuum spaces of the real world with the discrete spaces 
of the lattice world? 

The simple answer is that we take the number of lattice sites, N, to infinity. However, 
by itself, simply taking N — > oo does not provide a definite connection with the real world. 
Rather, we must also specify how the lattice spacing, a, is to behave as N — > oo. Now, as a 
general mathematical proposition, there is of course an infinite number of ways in which we 
can make a depend functionally on N. However, there are three particularly simple choices 
of a(N) that allow a direct physical interpretation, and are of general importance in all 
practical applications of Fourier theory. I will now discuss each of these in turn, and review 
explicitly how the finite-lattice results of Sec. 2 are transformed into their infinite-lattice 
counterparts. 

The simplest choice of a(N) is just 

a(N) = a = constant. (17) 

In other words, the real distance between lattice sites is kept constant; as we increase N, 
the extra sites simply increase the length of the lattice linearly, as shown by (11): 

L(N)=Na, i.e., L(N) oc N. (18) 

The constant spacing a remains a part of the mathematical description in the continuum 
limit N — > oo, and hence must correspond to some physical property of the system being 
modelled; in other words, a is a physically real quantity, rather than just a mathematical 
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symbol of our intermediate calculations. In practice, the prescription (17) would only be of 
use for systems that do, physically, have a lattice structure of periodicity a, such as a crystal. 

Eq. (13) tells us that the "length" of the momentum-space lattice will be a constant, 
regardless of the value of N: 

A(N) = — = constant. (19) 
a 

This implies that, as iV increases, more and more momentum-space lattice sites will be 
packed into this fixed "length" of A; Eq. (12) tells us that the "spacing" in momentum space 
between adjacent lattice sites is 

<») 

Hence, in the limit N — > oo, momentum space will become a continuum, within the "length" 
A = 2it/a. This is the "first Brillouin zone" of the crystallographer. 

We must now determine how the Fourier transform expressions of Sec. 2 are to be trans- 
formed in the limit N — > oo. Clearly, sums in position space will continue to be sums, but 
now over the infinite number of lattice sites. It is usual to first change the range of n from the 
1 to AT of Sec. 2 to a more "symmetrical" definition, such as from the least integer > —N/2 
to the greatest integer < +N/2, so that in the limit iV — > oo the sum over n will be from — oo 
to +oo. (The alternative would be to have an infinite lattice that possesses a single end, like 
a ray in geometry; this somewhat complicates the mathematical description, because instead 
of the continuous Fourier transform we would end up with the Laplace transform, which is 
not quite as convenient for our purposes.) 

On the other hand, sums in momentum space will clearly be converted into integrals, 
over the continuum of the first Brillouin zone. We have to be careful, however, in estab- 
lishing this correspondence. To illustrate why this requires extreme care, it is useful to first 
examine the definition of the convolution in 6-space, namely, Eq. (6). The usual transition 
to the continuum limit entails the insertion of a factor of 5p inside the sum, balanced by a 
compensating factor of 5p outside it: 



aVN 



(F*G) b = —-—'£8pF v G b - V = -^—^SpFvGb-v, (21) 
5pVN v 2tt v 

where it is understood that the sum over b' is over the range of integers of the first Brillouin 
zone. However, if we now take the limit N — > oo, we do not obtain the result we might 
expect: 



lim ^L^SpFvG^ = lim ^ [dp' F(p')G(p-p') = oo, (22) 

N^oo ZTX y iV— >oc Zn J 

where we are using (12) to relate btop (and b' to p') , and we are assuming that the integral of 
F(p')G(p—p') is finite. The problem is the factor of y/N outside the integral, which diverges 
as iV — > oo. 

To avoid this fate, we need to redefine all functions in momentum space in the limit 
iV — > oo, by multiplying them by a multiplicative factor that is a function of N. This is 
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essentially an elementary application of the idea of renormalisation. Now, since F(p) and 
G(p) appear in the convolution operation symmetrically, we must apply the same renormal- 
isation factor to each of them. However, the result of the convolution operation is itself a 
function in momentum space, and so should also be renormalised by the same factor. Taking 
these considerations into account, the divergence of (22) implies that, corresponding to any 
function H(jp) in momentum space, we should define a renormalised function 

H{p) = H(p)VN, (23) 

so that 



If we insert these renormalised functions into (21), we now obtain 



(F*G) b _ 1 y^ 5p F b ,G h 



so that 

1 



(F*G) b = — ]T>F b ,GU' = — ]T> F v Gw, 
We can now take the N — > oo limit without encountering any divergence problems: 

W^G)(P) = jJdp'F(p')G(p-p>), (25 ) 

A 

where the symbol U A" underneath the integral reminds us that we are integrating over a 
"box" of width A in momentum space, namely, the first Brillouin zone (regardless of how 
we choose its boundaries). The result (25) is the natural definition of the convolution of 
two continuous fields, as it corresponds to an "averaging" process, with the factor of 1/A 
balancing the integration over a box of width A. 

We can now investigate how the inverse Fourier transform (2) relates position space and 
momentum space in the case of constant a as N — > oo. Using the renormalised function (23) 
via Eq. (24), Eq. (2) gives 

f n = lim —L = J25 pe *«»n/Njb = = lim ± V<5pe 2 ™ 6 "^F 6 , (26) 



so that 



Thus, by using the renormalised function F(p) in momentum space, we have automatically 
obtained a sensible result for the inverse Fourier transform. Likewise, Eq. (1) gives 

F 1 - N/2 
4^ = -i= E e~ 2mbn/N f n , (28) 
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so that the limit N — > oo is similarly well-defined: 

oo 

F(p) = £ e_ ^ n /n- (29) 

n=— oo 

Finally, our use of renormalised functions in momentum space means that we need to refine 
the definition of the convolution operation in position space. Namely, if we form the product 
of two renormalised functions F(p) and G(p) in momentum space, then what is the corre- 
sponding operation between the functions /„ and g n in position space? The simplest way to 
obtain the correct result is to directly obtain the inverse Fourier transform of (FG)b- Using 
(26) without the introduction of 8p, together with (28), we obtain 



l _Y^ e ^ ibn / N F b G b = ^Y,e +2 ™^ /N Y, e ~ M ' N fn>Y, e ~ M ' /N 9n>> 

1 



N 



2mb(n-n'-n")/N 



— J^y^fn'^ 9n" N 6 n » n - n ' 
n' n" 

n' 

Thus, taking the limit N — > oo, and defining 

F-\{FG)(p)} = (f*g) n , 
we find that the position-space convolution operation is also completely well-defined: 

oo 

(f*g)n = E fn>9n-n>- (30) 
n'=— oo 

So far we have examined only the case in which the n — > oo limit has been taken with the 
lattice spacing kept physically constant, namely, Eq. (17), which gives us an infinite lattice 
in position space, and a continuum limited to the first Brillouin zone in momentum space. 
The second simple prescription we may make for the iV — > oo limit is to assume that the 
length of the position-space lattice, L, remains constant for all iV: 

L(N) — L — constant. (31) 

From (11), this implies that a(N) is inversely proportional to iV: 

a(N) = ± i.e., a(JV)«l (32) 

Eq. (13) now tells is that it is the "length" of the lattice in momentum space, A, that is 
proportional to N: 

A(N) = i.e., A(N) oc N, (33) 
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whereas it is the "spacing" in momentum space that is now a constant: 

Sp= 2 -^. (34) 

In other words, we have the complete conjugate of the previous scenario: we have a con- 
tinuum in position space, contained within a "box" of width L, which corresponds to an 
infinite lattice of states in momentum space. Because of the symmetry of the (finite lattice) 
Fourier transform, we automatically know that, in this case, we will have need to renormalise 
functions in position space: 

h(x) = h{x)VN, (35) 

so that 

h(x) = ^± (36) 



Now, using the fact that the spacing between lattice sites in position space is just 

Sx = a, 

the convolution operation in position space, (9), becomes 

1 



(f*9)n = J^J2 a fn'9n-n', 



n' 



which in the limit N — > oo gives 

(f*g)(x) = ^Jdxf(x')g{x-x'), (37) 

L 

where again the factor of 1/L naturally balances the integration over a box of width L; i.e., 
we are computing an average over the width of the box. 

Clearly, when we write the results in this way, the symmetry between position space and 
momentum space is manifest. This allows us to write down the constant-!/ analogues of (27), 
(29) and (30) directly: 

F b = ydxe~^ x f(x), (38) 

L 

oo 

f( x ) = ]T e-***F b , (39) 



-oo 
oo 



(F*G) b = J2 F b'G b ~ b >, (40) 



b'=— oo 

which can be confirmed by straightforward calculations. 
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Finally, we consider the case in which we want both position and momentum space to 
become continua in the limit N — ■> oo. Again, there are many ways in which this can be 
done, but the simplest to analyse is the symmetrical prescription 

a{N) K 7n' (41) 



so that the size of the "box" in both position space and momentum space is proportional to 
y/N: 

L(N) oc y/N, A(N) oc y/N, (42) 

and the spacing between adjacent lattice sites in both position and momentum space is 
inversely proportional to a/TV : 



6x(N) oc -L, Sp(N) oc -L. (43) 



The prescription (41) requires us to fix the constant of proportionality — say, 

a(N)= a ° 



N 

where ao is some constant with the dimensions of length. However, it is practically convenient 
for us to choose our system of units so that a = 1, so that (42) become 

L(N) = y/N, A{N) = 2nyfN, (42') 

and (43) become 

S X (N) = S p(N) = (43') 



If this is done, the iV — > oo results can be written down directly from (1), (2), (6) and (9): 

/oo 
dxe- ipx f(x), (44) 
-oo 

f( x ) = -L Tdpe^Fip), (45) 

Z7T J-oo 

1 r°° 

(F*G)(p) = — dp'F(p')G(p-p'), (46) 

ZTT J-oo 

[■OO 

(f*g)(x)= dx'f(x')g(x-x'), (47) 



with no need to "renormalise" functions in either position space or momentum space. It is 
an accident of pedagogical history — that should probably be rectified — that the continuum 
Fourier transform of Eqs. (44), (45), (46) and (47) is often taught to undergraduate students 
without providing a rigorous justification on the basis of it being the symmetrical iV — > oo 
limit (41) of the discrete Fourier transform. However, with the growing importance of lattice 
field theory, perhaps this prejudice will slowly be overcome. 
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For practical lattice calculations, of course, the dependence of a(N) on N does not fall 
neatly into any of the three simple categories listed above; indeed, one of the end-products 
of a lattice calculation is an estimate of a based on physically identifiable results. The 
a — > and L — > oo limits are still extrapolated, but they are not both tied symmetrically 
to N in the way described by (41); rather, both N and the gauge coupling j3 are varied 
from calculation to calculation. This does not at all invalidate the connection made above 
between the lattice and continuum descriptions, but merely represents a more complicated 
way of taking the double limit than the simple prescription (41). 

4. The problem for lattice field theory 

In Sec. 2 we looked at the relationship between a discrete lattice and its Fourier transform, 
from a purely mathematical point of view. In Sec. 3 we began to add some physical "meat" 
to the description, in terms of identifying position space and momentum space, but the 
results were still essentially mathematical in nature. In this section we will look at where 
the properties of Fourier theory become problematical for lattice field theory calculations in 
particular. 

The crux of the problem is the convolution operation. Even though we have set up 
our physical system to reside within a "box", in both position space and momentum space, 
the convolution operation goes outside this box and "reaches around the other side" of the 
lattice. Now, as noted above, this is not a problem in position space, because we are using 
the fiction of position space being wrapped into a ring (in higher dimensions, of spacetime 
being wrapped into an n-torus) in order to avoid having rigid boundaries that would destroy 
translation invariance and distort the calculations being performed inside the box. And, 
indeed, we find that for even the most elementary applications of a discrete space, we do 
wish to perform convolutions in position space; for example, the spatial derivative operator 
is really just a convolution, which corresponds to multiplying the function in question by ip 
in momentum space. 

It is with convolutions in momentum space that we have trouble. This corresponds 
to multiplying two functions in position space, which is something that is only done for 
nonlinear applications — like quantum field theory. Now, we must expect that conventional 
Fourier transform wisdom applicable for linear applications should need to be augmented 
when we delve into the nonlinear domain; and indeed this is the case. When we multiply two 
fields, the high-momentum components of each field combine to "alias" (or "masqerade") 
as low-momentum components in the result. A simple example will suffice to demonstrate 
the problem. Imagine we have a field in the highest possible momentum eigenstate in the 
positive-x direction, namely, p = +7i/a. Its spatial dependence will then be of the form e i7nE//a . 
Imagine that we now have another field, also in the highest possible momentum eigenstate 
in the positive-x direction. When we multiply the two fields in position space, the product 
then has a spatial dependence of simply e 2t7TX / a ^ which corresponds to a momentum value 
of +2n/a. Indeed, in general, when we multiply in position space two fields in momentum 
eigenstates of p\ and p 2 , the product is in a momentum eigenstate of p = p\ + p 2 - When one 
works it through, this simply represents the principle of conservation of momentum, when 
quantum field theory is interpreted in terms of particles of definite momentum. However, 
consider the lattice representation of these two fields. For spatial dependence 
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of e tnx/a simply becomes e lnn = (-1)" so that the field oscillates +1,-1, +1,-1,..., along 
the position-space lattice. If we multiply this by another field having the same spatial 
dependence, then we end up with +1,+1,+1,+1,.... In other words, the product of the 
two fields has vanishing momentum! Where did the momentum of 2ir/a disappear to? It 
was "absorbed by the crystal in an Umklapp process" — except that in quantum field theory 
we have no crystal. This zero- momentum eigenstate is a purely mathematical artefact. 

Another example might suffice to give a feel for the phenomenon of aliasing. Imagine that 
we have a field in a momentum eigenstate that is only 2/3 of the maximum momentum that 
can be represented on the lattice, namely, p = 2ir/3a, so that it has a spatial dependence of 
e 2i7rx/3a p or x ^ _ na ^ j g j us ^. e 2i7rn/3 means that its position-space dependence is 

1 1 1 1 

1, h 1 , I , 1, V % , I , .... 

2 2' 2 2' ' 22' 2 2' 

If we multiply this field by another field that is in the same momentum eigenstate of 
p = 2n/3a, then we again need to simply square each of these values: 

_1_.^3 1 V3 _I_-^3 1 , ->/3 

2*2' 2 +t 2 ' ' 2*2' 2 +l 2' ••" 

We now recognise this to be simply g-^WSa. ^ e . ; the product is in a momentum eigenstate 
of value p = —2ir/3a. In the language of quantum field theory, two particles of momentum 
p = +27r/3a have combined to create a particle of momentum p = — 27r/3a! 

Of course, the Umklapp analogy again tells us what is happening: the two particles com- 
bined to create a momentum eigenstate of p — +An/3a, and then the nonexistent "crystal" 
absorbed a momentum of 2n/a to bring the result back into the first Brillouin zone. 

Let us now consider how we might repair this oversight. 



5. A solution to the aliasing problem 

The engineer provides a simple methodology for solving the aliasing problem: ensure that 
the fields are filtered suitably, so that there is no "cross-talk" between aliased components 
and genuine components in momentum space. 

So how do we carry out such a programme? 

The obvious "knee-jerk" reaction would be to filter every field with a low-pass filter, 
allowing through only momentum components —n/2a < p < +n/2a that are no greater than 
half the maximum possible momentum (in absolute value) that can be represented on the 
lattice. When two such fields are multiplied together in position space, we can rest assured 
that the convolution process in momentum space will not "bleed out" of the first Brillouin 
zone; the sum of any two momenta can just reach the Brillouin zone boundary, but it cannot 
exceed it. This product of fields is then itself subjected to the same low-pass filtering, so that 
if we subsequently need to multiply it by some other (likewise filtered) field, the result again 
is guaranteed to be contained within the first Brillouin zone, without any aliasing artefacts. 

This process ensures that aliasing is removed completely, but it has cut down the available 
momentum modes too harshly — it is too conservative. Consider, instead, a low-pass filter 
that allows through all momenta — 2n/3a < p < +2n/3a that are no greater than two-thirds 
of the maximum possible momentum value (in absolute value) that can be represented on the 
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lattice. If we calculate the position-space product of two fields filtered in this way, we know 
that the result will contain momentum values that rightly belong as high up as p = +4n/3a 
and as low as p = — 4n/3a. The momentum components outside the first Brillouin zone will 
"bleed out" and be aliased in the first zone. Those momentum components that rightly 
belong in the range +ir/a < p < +47r/3a will be "aliased down" by the subtraction of a 
"crystal momentum" of 2n/a into the range — n/a < p < —2ir/3a. Likewise, those compo- 
nents that rightly belong in the range —in /3a < p < —n/a will be "aliased up" by addition 
of a "crystal momentum" of 2n/a into the range +2n/3a < p < +ir/a. In other words, the 
momentum ranges —ir/a < p < —2i\j3a and +27r/3a < p < +n/a are a mess: they are a 
combination of the genuine components of the product, polluted by the aliased components 
from the other end of the spectrum. 

This does not seem to be much of an advance. However, consider what we now do: we 
filter the result with the same low-pass filter, so that it is ready for any further products that 
we may need to form in position space. This filters out the polluted momentum ranges com- 
pletely; we are left with the components in the range —2ir/3a < p < +27r/3a that represent 
the product with perfect fidelity! (By "perfect fidelity" I do not mean that they represent 
the continuum exactly, because the lattice has provided a high-momentum regulator that 
removes momentum modes outside this range from the outset; rather, we have a perfect 
representation of the regulated theory, within this regulated momentum range.) 

In other words, the price we need to pay to remove aliasing from our lattice calculations is 
the complete elimination of one-third of the momentum modes in each dimension. For four- 
dimensional calculations, this means that we retain (2/3) 4 = 16/81 of the modes originally 
present on a lattice of any given dimensions, so that roughly four-fifths of them have been 
discarded (there are lots of "hyper corners" in a hypercube!). 

This reduction in momentum modes might seem to be wasteful: doesn't this mean that 
we will need to spend five times as long performing the same lattice calculation with the 
same hardware? I must emphasise that nothing could be farther from the truth. The 
momentum modes that we are eliminating are actually distorting calculations that are done 
without their elimination. Indeed, in recent years there has been a recognition of what has 
been called "noise" in the high-momentum components (short- distance structure) of lattice 
fields, and attempts have been made to "smooth" the fields slightly, to filter out some of 
this "noise". However, if this "noise" was truly noise — uncorrelated statistical hiss — then 
it would not affect lattice calculations systematically at all; it would simply add to the 
inherent noise of the Monte Carlo process, and disappear in the mean. Aliasing, in contrast, 
does not disappear on the average: it represents a form of "wrap-around" in momentum 
space, that can make finite contributions systematically. Indeed, since loop diagrams in field 
theory inherently diverge for high momentum, there is the potential for aliasing artefacts to 
completely swamp the genuine contributions to any calculation. 

6. Exact expressions for the antialiasing operator 

In the previous section I argued that fields in lattice field theory calculations should be 
passed regularly through an "antialiasing filter" to eliminate all aliasing from the resulting 
calculations. The required filter is what the engineers call a "low-pass filter": it allows 
through momentum components in the range — 2n/3a <p< +2n/3a without distortion, 



17 



and eliminates completely all components of higher momentum (in absolute value). 

The position-space form of such an operator for an infinite lattice can be derived by 
using the same method as was employed in [1] for the SLAC derivative operator, and the 
corresponding finite-lattice expressions can be derived by using the same method as was 
employed in [2] to "wrap" the operator around the finite lattice an infinite number of times. 
These are the tasks that we shall carry out in this section. 

The "trick" of [1] was to start with what we described in Sec. 3 as the "constant a" 
limit of iV — ■> oo, namely, we consider an infinite lattice in position space with spacing a, 
which corresponds to a continuum in momentum space within the first Brillouin zone. Our 
antialiasing is then simply described by 



F(p) 



1 if -27r/3a < p < +2?r/3a, 
otherwise. 



(48) 



The ( "renormalised" ) inverse Fourier transform (27) can be performed in closed form: 



Jn 



+27r/3a 



dpe 



tpx n 



^ e +2wix n /3a _ g - 



2nJ-2w/3a 2nix n 
which, from the definition of the sine function, and x n = na, yields 

a sin(27rn/3) 



roc 
Jn 



(49) 



(The superscript of "oo" reminds us that we are working with an infinite lattice at this stage.) 
Let us work in "lattice units" and set a = 1. To calculate fowe need to take a limiting process 
for n — > 0, for which sin(27m/3) rs 27rn/3, so that f = 2/3. For all other n = (mod 3), /„ 
is proportional to sin(27rA;) where k = n/3, and hence /„ vanishes. For n = 1 (mod 3), we 
have sin(27r/3 + 27r/c) = V / 3/2, so that f n = V3/2irn. Likewise, for n = 2 (mod 3), we have 
sin(47r/3 + 27r/c) = —V3/2, so that f n = —\Pij2-nn. Thus, putting these together, we have 



Jn 



2 
3 

V3 



2irn 
V3 



2im 




if n = 0, 

if n = 1 (mod 3), 

if n = 2 (mod 3), 

if n = (mod 3) and n^O. 



(50) 



We find that f n is an even function of n; it has a central term of 2/3, and as we move 
away from n = the terms alternate +, — , 0, +, — , 0, . . . , falling off as 1/n, with a numerical 
coefficient of V3/2tt. The sum of all the terms is 



^ fn 



2 V^vV 1 

3 7T ^ V3A; + 1 



3A; + 2 
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as expected: the filter f n spreads a delta function out across the lattice without changing 
its integrated value, and so when convolved with any other function g n will spread it out 
somewhat (filtering out the highest momentum states) without changing its average value. 

Let us now consider obtaining closed-form expressions for for a finite lattice of N 
sites. Since the general method has been described in lengthy detail in [1], I will not repeat 
it here, but merely cut directly to the intermediate results. 

Clearly, the simplest analysis is for the case when N = (mod 3), because the groups of 
three successive terms in the infinite-lattice result (50) will line up in columns. Specifically, 
all of the n = (mod 3) will fall on every third site, and so /„ will vanish for all n = (mod 3) 
except n = 0, which will retain its infinite-lattice value of 2/3. For n ^ (mod 3), we obtain 
the sum 

iV=0(mod3)_l . /2?m\~ J 1 1 



f iv=o(mod3) _ 1 (2irn\^ ( 

/„*0(mod3) ~ n Sm { 3 



Nk+n Nk + N-n 



\ 1 ■ 1 


'2im\ 




> = — sin 






J N 1 


v 3 j 





In actual fact, this expression is also correct for n = (mod 3), because for n^O the factor 
sin(27rn/3) vanishes, and the limit as n — > is the correct value of 2/3. Thus in full generality 
we may write 

fn ~» 4-(Hr) C0 O (51) 

where the limiting procedure for n = (yielding the value 2/3) is to be understood. 

For N^O (mod 3), the groups of three terms "roll through" every time we loop once 
around the lattice. The simplest way to proceed is to make use of this 3iV-periodicity rather 
than the iV-periodicity of the iV = (mod 3) case, and to sum up the two different sequences. 

Let us start with N = 1 (mod 3). The n = term picks up a correction: 

f MH3)J ^y(J 1 \ = l( l + J_\ (ro) 

Jn -° 3 Ti ^ {3Nk + N 3Nk + 2N) 3V 2NJ' v ; 

For n = 1 (mod 3) we obtain 

rN= 1 (mod 3) 

J«El(mod 3) r>_. / / 



(mods) " 2n j^\\3Nk+n 3Nk + 3N-n 



\3Nk+(N+n) 3Nk + 3N-(N+n) / 
1 f / nn \ / 7r nriX) 

= ^{ co \m- co \3 + m)\- 

For n = 2 (mod 3) we obtain 

rN= 1 (mod 3) 

Jn = 2 (mod3) 2tc ^ Q \\3Nk+(2N+n) 3Nk + 3N-(2N+n)J 

"( 1 



3Nk+n 3Nk+3N-n, 



2NV3 I 



f /2n 7in\ ( imW 
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Finally, for n = (mod 3), n ^ 0, we obtain 



f N= 1 (mod 3) 

Jn = o(mod3) 2tt |\3M;+(N+ ra) 3iVfc + 3iV-(iV+n) 



1 



^3iVA; + (2iV+n) 3iVA; + 3A^-(2iV+n) / 
1 f / 7r 7rn \ /27T ra N 

= 2^i C °H3 + 3wJ- COt lT + 3iV, 

We can now see that there are three sequences interleaved here: one with the cot of 7rn/3N, 
one with n/3 added to the angle, and one with 27r/3 added to the angle. We can thus 
combine them back into a single expression: 



, A ,-i , mA j 1 f . /2nn\ ( Tin \ . /2ir 2nn\ fix irn\ 

WV=l(mod3) = J in W +S i n 1_ COt - + 

Jn 3AM \ 3 J \3NJ \3 3 ) \3 3N J 



' 4tt 2irn\ f2ix im 



(53) 



where again we find that the n = term is automatically included as the n — > limit of the 
generic expression. 

Finally, we can carry through the same analysis for iV = 2 (mod 3). When all the terms 
are recombined, the only difference from (53) is that the "rolling" factor "rolls" in the reverse 
direction: 



1 f /27m\ / Tin \ / 4ir 2im\ ( tt Tin 

' 1 1 +sm 1 cot — | 

3NJ \ 3 3 J \3 3N 



/2vr 27rn\ f2n 7m\) ^ ^ 

+ sin 1 cot 1 )}. 

V 3 3 J V 3 3iVyi 



Again, we find that the n = term is included as the n — > limit of the generic expression, 
in agreement with the first-principles result: 



f JV=2(mod3) = 2 , V3y^f 1 1 \ = -( (tf) 

Jn=0 ~ 3 7r 1 3Nk + 2iV 3Nk + N) 3\ 2NJ' 



Finally, we can combine together the results (53) and (54) for N ^ (mod 3), by means 
of a simple factor in the argument of the "rolling factor" that depends on the value of 
k = N (mod 3): 



2kn 2ixn\ . /7r nn 
3 + 3N 



4&7T 27rn\ f2n im 



- sin j — - | cot { — 



3N 



(56) 



The expressions (51) and (56) are the final results of our calculation. These weights may 
be calculated for any value of the lattice size N, stored, and used to provide an antialiasing 
filter in that dimension. For a multi-dimensional lattice, we need simply perform the filtering 
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operation in each dimension separately. If the dimensions have different sizes, then we need 
to calculate and store a separate set of weights for each different size. Apart from that, 
however, the weights remain fixed for the entire run of the lattice calculation. 

7. Stochastic implementation of the antialiasing filter 

Clearly, the antialiasing filter described in the previous section is not suitable for direct 
implementation in practical lattice calculations. It is as horrifically nonlocal as the SLAC 
derivative operator studied in [1,2,3], linking a given site to every other site along the 
dimension of the lattice that we are filtering, with the absolute value of each term falling 
off only inversely with distance in the infinite-lattice case (with small numerical corrections 
for finite lattices). This means that it is not possible to truncate the series at some finite 
distance without severely distorting the fidelity of the operator (like we can with operators 
that fall off exponentially with distance) — which is precisely the fate that we want to avoid. 

Thus, any practical implementation of the antialiasing filter would almost certainly em- 
ploy the stochastic approach suggested in [1] for the SLAC operators. Namely, the absolute 
value of each weight calculated in the previous section will be taken to be the probability 
that the given linkage will be made between the site at which the antialiasing filter is being 
applied and the site which is at the given distance from the said site. Since the magnitude 
of the terms fall off inversely with distance, this implies that the average number of linkages 
that will be performed per site will be proportional to the logarithm of the length of the 
lattice in the given dimension, rather than being proportional to the length itself. 

This is precisely the philosophy espoused in [1, 2, 3], and so does not need to be elaborated 
on further here. However, I need to point out one important practical error contained in 
those previous papers, and one improvement in the algorithm that is of particular interest 
to the current case. 

The error relates to the requirement that each and every step of a lattice calculation 
maintain the unitarity of the formalism. In each of [1,2,3] I either stated or implied that 
unitarity would be maintained if every computation of a positive-distance term from a given 
site was always accompanied by the computation of the corresponding negative- distance term 
from that given site. This erroneous claim was based on a misapplication of the requirements 
of Hermiticity of the momentum operator. When the Hermiticity requirements are examined 
in closer detail, what we find is something that looks somewhat similar at a glance, but 
which has completely different ramifications for practical calculations, namely, that if for a 
Hermitian operator we perform a calculation at any site i that depends on the value of a 
field at any other site j, then we must perform the Hermitian-conjugate calculation at site j 
that depends on the Hermitian-conjugate value of that same field at site i. In other words, 
we have a form of "reciprocity requirement" : if a site "gives", then it must "receive". It is not 
at all required that a positive-distance calculation from a given site be always accompanied 
by a negative-distance calculation from that same site (although this will be true on the 
average). 

The algorithmic improvement relates to the implementation of the probabilistic aspect of 
each possible "linkage" between sites. In [1] I pointed out that the evaluation of a probability 
of 1/n (with n an integer) can be carried out much more efficiently than for an arbitrary 
real probability between and 1, because the former can be determined directly from the 
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pseudo-random bitstream, without the need for floating-point calculations at all. This was 
very fortuitous for the considerations of [1], because the absolute value of each term of the 
first-derivative operator for an infinite lattice is precisely 1/n, where n is the distance in 
lattice sites. This aspect of my argument remains true. 

However, when in [2] I "corrected" the infinite-lattice results for finite lattices, the ab- 
solute value for the coefficient for a distance of n was changed from 1/n to some other 
number that, for small n, was very close to 1/n. I therefore proposed that the probability 
for performing each calculation be maintained at 1/n, because this yes/no decision needs to 
be performed as quickly as possible, as it performed for every distance from the given site. 
Rather, if the result of this probabilistic decision was to actually perform the linkage, then 
the calculation would be multiplied by a factor that was the ratio of the true (corrected) 
weight to 1/n. The argument goes that, since these "correction factors" are all close to 
unity, the optimal balance obtained in [1] between the statistical noise of the operator and 
the minimisation of the number of calculations is effectively maintained, without the need 
to introduce floating-point calculations into the innermost probabilistic decision loop. 

The logic of this argument also remains true. However, what I overlooked at the time 
was the fact that, as we approach the maximum possible distance from the site in question, 
the "correction factors" can vary substantially from unity. Of greatest practical importance 
is the case of the first-derivative operator for lattices with an even number of sites, because 
this is the situation we are usually considering for lattice field theory. In this case, the 
"correction factors" actually vanish as we approach the maximum distance. This means 
that, in an "ideal" implementation (which wouldn't care about the cost of floating-point 
operations), the probability of performing a calculation at large n will fall off much faster 
than the 1/n of the infinite-lattice operator. By maintaining the probability at 1/n, we 
would be computing a greater number of linkages than we need to: the extra linkages would 
come in with a weight substantially less than unity, which means that they are wasteful of 
computing power when we consider the balance between statistical noise and the number of 
computations performed. 

The solution is to slightly change the philosophy of the computation of the probabilistic 
weight and the correction factor. Namely, we forget about the particular structure of the 
SLAC derivative operator altogether, and consider the case of any arbitrary probability p, 
where < p < 1. We now compute M = round(l/p), where round(rr) is the closest integer 
to x. In the innermost lattice loop, we efficiently generate a random integer between 
and M — 1, as described in [1]. If this random integer vanishes, we perform the linkage 
calculation, and multiply it by a weight of Mp. This process then ensures that we are always 
performing a linkage that has a weight that is close to unity, and hence is maximally efficient. 
(The greatest deviation from unit weight occurs for high probabilities: For p > 2/3, we find 
that M = 1 and hence the linkage is always performed, with a weight that is less than 
unity, but no less than 2/3. For 2/5 < p < 2/3, we find that M = 2, so that the calculation 
is performed half the time, with a weight that can lie in the range 4/5 < Mp < 4/3. For 
p < 2/5 the weight will clearly not vary from unity by more than 20%.) 

Simulations of the stochastic version of the SLAC derivative operator, using this ap- 
proach, have successfully demonstrated the expected behaviour. The Fourier transform of 
the operator has a mean value of precisely ip, with statistical noise that is distributed almost 
completely uniformly between all of the momentum-space components (except p = 0, which 
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vanishes by symmetry for odd N, and has a slightly larger value than the other momentum 
components for even N). 

As a final note, the modified probabilistic algorithm described above allows us to "dial 
the stochasticity" of the given operator, if desired. To do this, we modify the computation of 
M to M= round(r] / p) . If the resulting M vanishes, then we set it equal to unity. If we set 
the dial to rj — 1, we obtain the probabilistic algorithm described above. If we set the dial 
to i] = 0, then we will always obtain M — 1 (pipped up from M — 0), and so we will always 
perform the linkage, which means that we have removed completely the stochastic nature of 
the operator. For < rj < 1 we have an intermediate situation, in which the statistical noise 
is less than for r) — 1, but not zero. Thus the parameter 77 is effectively the "stochasticity", a 
number that we can dial. We can even set rj > 1, which will perform linkages even less often, 
at the expense of correspondingly increasing the weights above unity, so that the statistical 
noise is increased. 
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